#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gmp.h>
#include <time.h>

typedef struct{
    mpf_t a;
    mpf_t b;
    mpf_t p;
}Data;

int main(){
    time_t start, end;
    Data D[23];
    int i;

//    start = time(NULL);
    mpf_set_default_prec(40000016);
    mpf_t v_sqrt;
    mpf_t temp;

    for(i=0; i<23; i++){
        mpf_init_set_ui(D[i].a, 0L);
        mpf_init_set_ui(D[i].b, 0L);
        mpf_init_set_ui(D[i].p, 0L);
    }

    mpf_init_set_ui(v_sqrt, 0l);
    mpf_init_set_ui(temp, 0L);

    mpf_sqrt_ui(v_sqrt, 2L);

    mpf_add(D[0].a, D[0].a, v_sqrt);
    mpf_add_ui(D[0].p, v_sqrt, 2L);

    for(i=0; i<22; i++){
        mpf_sqrt(v_sqrt, D[i].a);
        //Update of An+1
        mpf_ui_div(temp, 1L, v_sqrt); // temp = 1 / sqrt(A(n))
        mpf_add(D[i+1].a, v_sqrt, temp); // A(n+1) = sqrt(A(n)) + (1/sqrt( A(n)))
        mpf_div_ui(D[i+1].a, D[i+1].a, 2L); // A(n+1) = (sqrt(A(n)) + (1/sqrt(A(n)))) / 2


        //Update of Bn+1
        mpf_add_ui(temp, D[i].b, 1L); //temp = 1 + b(n)
        mpf_mul(D[i+1].b, temp, v_sqrt); //b(n+1) = (1 + b(n)) * sqrt(A(n))
        mpf_add(temp, D[i].a, D[i].b); //temp = a(n) + b(n)
        mpf_div(D[i+1].b, D[i+1].b, temp); // B(n+1) = (1 + b(n)) * sqrt(A(n)) / a(n) + b(n)


        //Update of Pn+1
        mpf_add_ui(temp, D[i+1].a, 1L); //temp = A(n+1) + 1
        mpf_mul(temp, temp, D[i].p); // temp = (A(n+1) + 1) * p(n)
        mpf_mul(D[i+1].p, temp, D[i+1].b); // P(n+1) = ( (A(n+1) + 1) * p(n) ) * b(n+1)
        mpf_add_ui(temp, D[i+1].b, 1L); // temp = (1 + B(n+1) )
        mpf_div(D[i+1].p, D[i+1].p, temp); // P(n+1) = ( (A(n+1) + 1) * p(n) ) * b(n+1) / (1 + B(n+1) )

        //Clear useless variables
        mpf_clear(D[i].a);
        mpf_clear(D[i].b);

    }

    mpf_out_str(stdout, 10, 40000016, D[i].p);
    
	return 0;
}
